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1 Introduction 

The recognition that the quark-gluon plasma (QGP) produced in relativistic heavy ion 
collisions is strongly coupled [1] has prompted much work using gauge/gravity duality 
(or “holography”) to study aspects of non-equilibrium dynamics in strongly coupled 
A/" = 4 supersymmetric Yang-Mills theory (SYM), which may be viewed as a theoreti¬ 
cally tractable toy model for real QGP. Work to date has involved various idealizations 
(boost invariance, planar shocks, imploding shells) which reduce the dimensionality of 
the computational problem, at the cost of eliminating significant aspects of the real 
collision dynamics [2-7]. In this paper, we report the results of the first calculation of 
this type which does not impose unrealistic dimensionality reducing restrictions. We 
study the collision, at non-zero impact parameter, of bounded, localized distributions 
of energy density which mimic colliding relativistic Lorentz-contracted nuclei and form 
stable incoming projectiles in strongly coupled SYM. 

In the dual gravitational description our initial state consists of two localized in¬ 
coming gravitational waves, in asymptotically anti-de Sitter (AdS) spacetime, which 
are arranged to collide at a non-zero impact parameter. The precollision geometry 
contains a trapped surface and the collision results in the formation of a black brane. 
We numerically solve the full 5D Einstein equations for the geometry during and after 
the collision and report on the evolution of the SYM stress-energy tensor . 
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2 Gravitational formulation 


2.1 Single shocks 

We construct initial data for Einstein’s equations by combining the metrics describ¬ 
ing gravitational shock waves moving at the speed of light in opposite directions. In 
Fefferman-Graham (FG) coordinates, the metric of a single shock moving in the ± 2 : 
direction is given by [8, 9] 

ds^ = — [—dt^ + dx\ -\- dz^ -\- ds^ + h±{x±, z=p, s) dz"^] , (2.1) 

s 

(where x± = {x, y} and 2 =p = z^t). This is a sourceless solution to Einstein’s equations 
with cosmological constant A = provided the function h± satisfies the linear 

differential equation 

(a2-fa3 + Vi)h± = 0. (2.2) 

Solutions to this equation which vanish at the boundary, s = 0, may be written in the 
form 

/ r/^k — 

g*kx.y± H±{k±, %) S{s‘^/kl) hik^s ), (2.3) 

where H± is an arbitrary function of the 2D transverse wavevector k± and the longi¬ 
tudinal variable 2 =^ (and I 2 is a modified Bessel function). 

The geometry described by eqs. (2.1) and (2.3) represents a state in the dual SYM 
theory with a stress-energy tensor expectation value given by [10]. 

(TOO) = (T^^) = ±(T°^) = kH^{x^,z^), (2.4) 

with all other components vanishing. Here, H± is the 2D transverse Fourier transform 
of i7±, and the constant k, = /{4:7:Gn) — N^/{2n^) (with Nc the gauge group rank 

of the SU{Nc) SYM theory). In other words, the function H± specifies the energy 
density (and longitudinal stress and momentum density) of a shock wave moving, non- 
dispersively, at the speed of light in the ±z direction. Given any choice of this energy 
density profile, eqs. (2.1) and (2.3) give an explicit form for the unique dual gravitational 
geometry which describes this shock wave. 

For simplicity, we consider Gaussian energy density profiles, 

H±{x±, z^) = exp {-Izl/w^) exp [-lix± T b/2f/R^] , (2.5) 

V 27rtc^ 

with longitudinal width w, transverse width R, and a transverse offset ±b/2. Hence, b 
will be the impact parameter when these oppositely directed shock waves collide. The 
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amplitude A gives the maximum value of the longitudinally integrated energy density 
(divided by k). 

For the explicit computations presented below, we adopt units in which the ampli¬ 
tude A equals unity and choose longitudinal and transverse widths re = | and i? = 4, 
respectively, and impact parameter b = ^Rx. Consequently, the stress tensor (2.4) 
describes localized lumps of energy centered about x = :^bf2, y = 0, and 2 ; = ±t. The 
AdS curvature scale L is not a physical scale in the dual QFT and may independently 
be set to unity. 

2.2 Infalling coordinates 

Our time evolution scheme [ 11 ] for asymptotically-AdS gravitational dynamics uses 
infalling Eddington-Finkelstein (EF) coordinates in which the spacetime metric has 
the form 

^2 

ds^ =—gnv{x,r)dx^dx'' + 2dr dt ^ (2.6) 

with Greek indices denoting spacetime boundary coordinates, = {t,x,y,z), and 
r the bulk radial coordinate. For the asymptotically-AdS geometries of interest, the 
metric coelRcients have the near-boundary behavior = 77 ^^ -|- gl^}/A + 0{l/r^) 

as r ^ 00 , with = diag(—1, -|-1, -|-1, -|-1) the usual Minkowski metric tensor. The 

sub-leading coefficients gl^J determine the SYM stress tensor expectation value. In 
these coordinates, 

{T^^)/k = gf} + \ glS- (2.7) 

Although the single shock solution to Einstein’s equations has the nice analytic 
form (2.4) in Fefferman-Graham coordinates, this geometry does not have a simple 
analytic form in infalling EE coordinates; the transformation to infalling coordinates 
must be performed numerically. One may easily show that, in infalling coordinates, 
curves along which r varies, with the other coordinates held fixed, are infalling null 
geodesics (with r an affine parameter). To transform the geometry (2.1) to the infalling 
form (2.6) one must locate the same congruence of infalling radial null geodesics in 
FG coordinates. Let Y = { 7 /^, s} denote the FG coordinates of some event, and let 
X = {x^,r} denote the EF coordinates of the event at affine parameter r along the 
radial infalling geodesic which begins at boundary coordinates x^. Then the solution 
Y(X) to the geodesic equation (in FG coordinates) for the same null geodesic which 
begins at boundary coordinates x^ provides the required mapping between EF and 
FG coordinates. Given this mapping, the required transformation of the FG metric 
components Gmn{^) to the metric components Gmn{^) in our infalling coordinates 
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(with hs2 = GMNiX) dX^dX^ = Gmn{Y) dY^dY^) is simply 

QyA QyB _ 

Gmn(A') = —— Gab(Y(X)) . (2.8) 

To compute the congruence Y{X), we first periodically compactified spatial direc¬ 
tions (to obtain a finite computational domain) with transverse size = Ly = 32 
and longitudinal length = 12. We employed spectral methods [11], and used a 
rectangular grid built from single domain Fourier grids with 32 points in the transverse 
directions, a 256 point Fourier grid in the longitudinal direction, and three Chebyshev 
domains of 32 points each in the radial direction. We used a Newton iterative procedure 
to solve the non-linear geodesic equation starting at each grid point on the boundary 
(modulo the C' 4 „ transverse cubic symmetry).^ We integrated the geodesic equations 
to a depth of s = 5 and fixed the residual radial shift reparameterization freedom [11] 
by demanding that the surfaces u = l/r = 5in infalling coordinates, and s = 5 in FG 
coordinates, coincide. 

2.3 Colliding shocks: initial data 

For early times, t <C —w, the Gaussian profiles H± of the oppositely directed incoming 
shocks have negligible overlap and the precollision geometry can be constructed by 
replacing the last term in the single shock metric (2.1) with the sum of corresponding 
terms from left and right moving shocks. The resulting metric satisfies Einstein’s 
equations, at early times, up to exponentially small errors. 

Initial data for the subsequent time evolution consists of the values of the spatial 
metric coefficients, rescaled to have unit determinant, gij = gij/{det ||fifq||)^^^, on every 
point of the initial time slice, together with the values of the energy and momentum 
density (or equivalently the asymptotic coefficients g^^J) at the initial time. The initial 
time slice was chosen to lie at t = —2, We slightly modified the initial data by adding 
a small uniform background energy density, equal to 3.7% of the peak energy density 
of the incoming shocks. This modestly displaced the location of the apparent horizon 
toward the boundary, improving numerical stability and allowing use of a coarser grid, 
thereby reducing memory requirements.^ 

^To obtain solutions with high accuracy, we used 40 digit arithmetic in this step. This allowed us 
to reduce both numerical arithmetic and iterative convergence errors to negligible levels, leaving the 
spectral truncation error as the limiting numerical issue. 

^Note added in proof: we have been able to recompute the evolution using a much smaller back¬ 
ground energy density, equal to 0.05% of the peak energy density. Over the time duration studied, the 
impact on the resulting dynamics of this change in background energy density is negligible. None of 
the features of the results discussed below are significantly affected by this background energy density. 
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2.4 Colliding shocks: time evolution 

Time development of the geometry was calculated using the procedure described in 
detail in ref. [11]. Time evolution was performed using a spectral grid of size = 
Ny = 39, Nz = 145, and Nr = 40. A fourth-order Runge-Kutta time integration 
algorithm was used with a time step of 0.005. Integration continued to a final time of 
t = 4. Even though our gravitational dynamics is five dimensional, with no symmetry 
restrictions imposed, we were able to perform the time evolution on a 6 core desktop 
computer with a 14 day runtime. 

3 Results 

In Fig. 1 we plot the energy density (top) and energy flux (bottom) in the 
plane y = 0 at several values of time.^ Note that the color scaling varies from plot to 
plot. The color scaling in the lower plots denotes |r°*| and the flow lines indicate the 
direction of At time t = —2 the system consists of two well separated lumps of 
energy moving towards each other at the speed of light. The non-zero impact parameter 
in the a:-direction is apparent. At time t = 0, when both shocks are centered at .2 = 0, 
the energy density and flux differ by only 8% and 15%, respectively, from a linear 
superposition of two unmodified shocks. Nevertheless, the subsequent evolution is very 
different from two unmodified outgoing shocks: the remnants of the initial shocks, 
which remain close to the lightcone, 2 ; = ±t, are significantly attenuated in amplitude 
with the extracted energy deposited in the interior region. 

At sufficiently late times, in accord with fluid/gravity duality [12, 13], the evolution 
of the stress tensor should be governed by hydrodynamics. To compare to hydrody¬ 
namics, we define the fluid velocity to be the normalized time-like = “1) future 

directed eigenvector of the stress tensor, 

(3.1) 

with e the proper energy density. Given the flow field and energy density e, we 
construct the hydrodynamic approximation to the stress tensor using the consti¬ 

tutive relations of first order viscous hydrodynamics. 

In Fig. 2 we plot the stress tensor components T** and and their hydrodynamic 
approximations Thy^ro ^hydro’ spatial origin x = y = z = 0, asa function 

of time. At this point the stress tensor is diagonal, the flow velocity w = 0, and 
and are simply the pressures in the x and 2 ; directions. As shown in the figure, the 

^Here and henceforth, is really the expectation value of the SYM stress-energy tensor divided 
by « = N2/(27r2). 
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t = -2 


t = 0 


t = 2 


t = 4 



Figure 1. The energy density (top) and energy flux |T^^| (bottom), at four different 
times, in the plane ^ = 0. Streamlines in the lower plots denote the direction of energy flux. 
Note that the color scaling varies from plot to plot. At the initial time t = —2 the shocks 
are at z = ±2. The non-zero impact parameter in the x-direction is apparent. The shocks 
move in the ±z direction at the speed of light and collide at t = z = 0. After the collision the 
remnants of the initial shocks, which remain close to the lightcone, z = ±t, are significantly 
attenuated in amplitude with the extracted energy deposited in the interior region. The 
development of transverse flow is apparent at positive times. 


pressures increase dramatically during the collision, reflecting a system which is highly 
anisotropic and far from equilibrium. After a time t 1.25 the pressures are well 
described by viscous hydrodynamics. Remarkably, at this time the transverse pressure 
Txx -g ten times larger than (This latter phenomena has also been seen in 

1 + 1 dimensional flow [2-4].) 
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Figure 2. Stress tensor components and at the spatial origin, x = y = z = 0, as 
a function of time. Dashed lines denote the viscous hydrodynamic approximation. Around 
t = 0 the system is highly anisotropic and far from equilibrium. Nevertheless, at this point 
in space, the system begins to evolve hydrodynamically at t ~ 1.25. 

To quantify the domain in which hydrodynamics is applicable, we define a residual 
measure 

A = (l/p)v/AT^,ArM-, (3.2) 

with AT^’^ = and p = ej3 the average pressure in the local rest frame. 

The quantity A is frame-independent but, when evaluated in the local fluid rest frame, 
reduces to the relative difference between the spatial stress in and Regions 

with A <C 1 are evolving hydrodynamically. 

In Fig. 3 we plot A in the transverse plane at proper times r = 1, 1.25, and 2, 
and rapidities ^ = 0 and 1. The color scaling is the same in all plots. Focusing first on 
^ = 0 (top row), at r = 1 one sees that A > 0.5 in the central region {x,y ~ 0), and 
hydrodynamics is not a good description. However, by r = 1.25 a fluid droplet with A < 
0.15 and transverse radius x± = |a;x| < 5.3 has formed, with subsequent evolution well 
described by hydrodynamics. At r = 2 the transverse size of the droplet has increased 
and A < 0.15 for x± < 8.6. Turning now to the behavior at rapidity ^ = 1 (bottom 
row), one sees that for small x± the system is already evolving hydrodynamically at 
T = 1. Moreover, the onset of hydrodynamics occurs earlier for a: < 0 than for x > 0. 
This feature reflects the fact that the receding maxima remain far from equilibrium 
and non-hydrodynamic, and (as seen in Fig. 1), the maxima with ^ > 0 lies at x > 0. 
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Figure 3. The residual A in the transverse plane, at several proper times r and two values of 
rapidity Regions with A <C 1 are evolving hydrodynamically. At ^ = 0 (top row) the central 
region becomes hydrodynamic at r 1.25, whereas at ^ = 1 (bottom row) hydrodynamic 
behavior of the central region has already begun by r « 1. At ^ = 1, hydrodynamic behavior 
first sets in at a: < 0. This feature reflects the fact that the receding maxima remain far from 
equilibrium and non-hydrodynamic, and the maxima with ^ > 0 lies at a; > 0. 

Interestingly, the inclusion of transverse dynamics seems to hasten the approach to 
local equilibrium: the equilibration time thydro 1-25 is about 30% smaller than was the 
case in our previous studies [3, 11] of planar shock collisions. Recent work [14, 15] has 
found that equilibration time scales of far-from-equilibrium states can be understood, 
at least semi-quantitatively, in terms of the spectrum of quasinormal modes. Post¬ 
collision, a distribution of quasinormal modes will be excited. With the inclusion 
of transverse dynamics, the dominantly excited quasinormal modes will be ones with 
non-zero transverse wavevectors of order 1/R. Faster decay rates of quasinormal modes 
with increasing |A:u| should lead to faster apparent relaxation of the entire sum (i.e., 
the metric perturbation) independent of location in the transverse plane. This will be 
the case provided one considers sufficiently coarse measures of relaxation, such as our 
A < 0.15 criterion, which can become satisfied when numerous quasinormal modes 
are comparably excited. With a much more stringent local equilibration criterion, we 
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z—x plane 


z—y plane 


x—y plane 



Figure 4. The fluid 3-velocity |'i;| at time t = 4, in the z—x, z—y, and x—y planes. Stream¬ 
lines denote the direction of v. Regions in which the residual A, deflned in Eq. (3.2), 
is greater than 0.15 have been excised; within these regions the system is behaving non- 
hydrodynamically. The maximum of |u|, which occurs in the vicinity of the receding maxima, 
is 0.64. In contrast, the maximum transverse velocity in the x—y plane is 0.3. 


would expect to see less of a difference relative to the planar case. 

A striking feature of the the post-collision evolution in Fig. 1 is the appearance of 
flow in the transverse plane at early times. The early-time acceleration imparted on the 
transverse flow can have a significant impact on the subsequent transverse expansion. 
In Fig. 4 we plot the fluid 3-velocity v = u/u^ in the z—x, z—y, and x—y planes at time 
t = 4. The color scaling, which indicates |t;|, is the same in each plot. The flow lines 
show the direction of v. Regions in which A > 0.15, and the system is not behaving 
hydrodynamically, have been excised. Already at time t = 4 and radius ~ 5 the 
transverse fluid velocity in the x—y plane has magnitude 0.3. In contrast, the maximum 
of the longitudinal velocity, which occurs in the neighborhood of the receding maxima, 
is 0.64. 

One sees from Fig. 4 that the fluid velocity in the x—y plane is nearly radial: we 
see no strong signatures of elliptic flow. At vanishing rapidity (or z = 0), the traceless 
part of the spatial stress, as well as the hydrodynamic residual shown in the upper row 
of fig. 3, deviate significantly from rotational symmetry in the x—y plane (at a 5-6% 
level at t = 4). However, the energy density and the resulting fluid velocity deviate only 
slightly from azimuthal symmetry; at t = 4 their transverse plane anisotropy is about 
1%. Our computed evolution does not extend through the entire hydrodynamic phase 
of the plasma and continued hydrodynamic expansion may generate more elliptic flow. 
However, the very modest elliptic flow, despite the substantial impact parameter, may 
reflect an unphysical aspect of our choice of Gaussian initial energy density profiles. 
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t = 1.25 t = 2 



Figure 5. The average lab-frame radial energy flux, as a function of x±, at two values each of 
proper time and rapidity. Also shown is the approximate form (3.4) which, for small rapidity, 
agrees quite well with the full results. 

Gaussian profiles differ, of course, from more realistic models of nuclear density (such as 
Wood-Saxon profiles). But a Gaussian profile in the transverse plane is computationally 
convenient. The very rapid decrease of its Fourier transform allowed us to use a rather 
modest transverse grid. However, a special feature of Gaussian transverse plane profiles 
is that the overlap function (i.e., the product of the energy densities of the two incident 
projectiles) is rotationally symmetric for any impact parameter: 

g-|Gx+b/2)/ii2 g-|(a;x-6/2)/iJ2 ^ ^_(^2^+(6/2)2)/ij2 ^ 

With Gaussian profiles (and equal size incident projectiles), there simply is no “almond” 
in the overlap function. The actual dynamics depends, of course, on the full spacetime 
dependence of the initial data, not just on this overlap function and, as noted above, 
the stress tensor and the residual A depart significantly from transverse rotational 
symmetry. Nevertheless, a simple picture in which the early transverse flow reflects, in 
large part, this initial overlap function appears plausible. 

A notable feature in our results is that the fluid flow, in the z—x plane, is not 
symmetric about the 2 : axis and the longitudinal flow does not vanish at 2 = 0. The 
latter observation is a direct violation of the simplified model of boost invariant flow, in 
which = zjt and vanishes at z = 0. Nevertheless, at t = 4 and in the region of space 
where e > 0.6max(e), the longitudinal flow is roughly described by boost invariant 
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flow at the 20% level or better, with larger deviations appearing at larger rapidities^ 
For planar shock collisions, the deviation of the longitudinal fluid velocity from boost 
invariant flow decreases as the shock thickness decreases [4, 11].^ It will be interesting 
to see if this also holds when transverse dynamics is included. 

We conclude by discussing the early-time transverse flow predicted in ref. [16]. 
There, using assumptions of boost invariance and transverse plane rotational symmetry, 
it is argued that at early times the transverse energy flux is proportional to the gradient 
of the energy density and grows linearly with time, 

= T^y = -ldye. (3.4) 

In Fig. 5 we plot the angular averaged radial flow together with the 

approximation (3.4), at proper times r = 1.25 and 2, and rapidities ^ = 0 and 1. The 
approximation (3.4) works remarkably well at both times and rapidities, although the 
agreement is not quite as good at ^ = 1 where the assumption of boost invariance 
is more strongly violated. It would be interesting to see if the agreement with the 
approximation (3.4) improves when the shock thickness decreases. 

4 Final remarks 

We have presented results from the first calculation, using numerical holography, of the 
evolution in strongly coupled A/" = 4 SYM theory of an initial state which resembles two 
colliding heavy nuclei, Correctly treating the dynamics, both longitudinal and trans¬ 
verse, requires solving a five dimensional gravitational initial value problem. This was 
challenging, but feasible, using our characteristic formulation and (good) desktop-scale 
computing resources. The results presented above are a proof-of-principle. Clearly, it 
will be interesting to explore how results change as the impact parameter, longitudinal 
thickness, and transverse size of each projectile are varied. It should also be possible 
to explore the effects produced by using more realistic initial energy density profiles, 
possibly including shorter distance fluctuations of the type which have been suggested 
to be relevant in real heavy ion collisions. (See, for example, ref. [17-20] and references 
therein.) We hope future work will shed light on some of these topics. 

^Specifically, in the region where e > 0.6max (e), we find max — z/t\ < 0.20 max |n^|. 

^However, even in the limit of thin planar shocks the proper energy density has strong rapidity 
dependence [4, 11]. 
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